#### Daniel J. Gustafson
#### 10 January 2017
#### gustafson@unc.edu

rm(list=ls(all=TRUE))

## Load packages
library(ggplot2)
library(xtable)
library(stargazer)
library(Zelig)
library(multiwayvcov)
library(lmtest)
library(plyr)
library(zoo)
library(reshape)
library(pastecs)
library(texreg)
library(scales)
library(car)
library(MASS)

## This is a truncated version of the MAROB data containing only variables for this analysis.
MAROB.Negot <- read.csv("/Users/danielgustafson/Documents/Grad/MAROB Negotiation/MAROB_Negot_new.csv")

MAROB.Negot <- na.omit(MAROB.Negot)

## Attach the data to call it by name
attach(MAROB.Negot)

## This is our model.
logit.1 <- glm(ACTIVENEGOTIATION ~ ORGST9 + orgmilitant + domorgprot + orgpop + relorg + natorgnew + polity2 + ORGST9*orgmilitant, data = MAROB.Negot, family = "binomial")

## Clustered standard errors
cluster.vcv <- cluster.vcov(logit.1, cluster = orgId)

## Results
new.coef <- coeftest(logit.1, cluster.vcv)

summary(logit.1)

vif(logit.1)
######################################
## This section contains code for interpreting the interaction term.

## Name your vcv. If you haven't done an clustering, the base command is vcov(model).
VCV.logit.1 <- cluster.vcv

## Number of coefs
P <- length(coef(logit.1))

## Set your CI
CI <- .95

## Values of militancy. Number of values = N, here.
s.orgmilitant <- seq(min(orgmilitant),max(orgmilitant),length.out=1135)

## The marginal effect of territorial control. 2 is the coef on terr. cont. 9 is the interaction coef. N.B. interactions are always listed last in coef lists.
Marg.ORGST9 <- coef(logit.1)[2] + coef(logit.1)[9]*s.orgmilitant

## Calculate the standard error for terr. cont.
SE.Marg.ORGST9 <- sqrt(VCV.logit.1[2,2] + (s.orgmilitant^2)*VCV.logit.1[9,9] + 2*s.orgmilitant*VCV.logit.1[2,9])

## Compute confidence intervals
Lowerci.ORGST9 <- Marg.ORGST9 - qt((1-CI)/2,1135-P)*SE.Marg.ORGST9
Upperci.ORGST9 <- Marg.ORGST9 + qt((1-CI)/2,1135-P)*SE.Marg.ORGST9

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

ggplot() +
labs(x="Observed Values of Militancy", y="Marginal Effect of Territorial Control") +
geom_line(aes(s.orgmilitant,Marg.ORGST9, color="Marginal Effect")) +
geom_ribbon(aes(x=s.orgmilitant, ymin=Lowerci.ORGST9, ymax=Upperci.ORGST9, fill="95% Confidence Interval"), alpha=.5) +
theme_bw() +
theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),legend.position="bottom") +
scale_fill_manual(values="black",name="")+
scale_color_manual(values="black",name="")

######################################

s.ORGST9 <- seq(min(ORGST9),max(ORGST9),length.out=1135)

## The marginal effect of militancy for varying values of territorial control.
Marg.orgmilitant <- coef(logit.1)[3] + coef(logit.1)[9]*s.ORGST9

SE.Marg.orgmilitant <- sqrt(VCV.logit.1[3,3] + (s.ORGST9^2)*VCV.logit.1[9,9] + 2*s.ORGST9*VCV.logit.1[3,9])

Lowerci.orgmilitant <- Marg.orgmilitant - qt((1-CI)/2,1135-P)*SE.Marg.orgmilitant
Upperci.orgmilitant <- Marg.orgmilitant + qt((1-CI)/2,1135-P)*SE.Marg.orgmilitant

## Plot for visualization. N.B. the only points that exist on this plot occur at territorial control = 0,1,2. This figure will not be presented.
plot(s.ORGST9,Marg.orgmilitant,type="n",
     xlab="Observed Values of Territorial Control",
     ylab="Marginal Effect of Militancy",
     ylim=c(min(Lowerci.orgmilitant)-2,max(Upperci.orgmilitant)+2))

lines(s.ORGST9,Marg.orgmilitant,lwd=5,lty=1)
lines(s.ORGST9,Lowerci.orgmilitant,lwd=4,lty=2, col="blue")
lines(s.ORGST9,Upperci.orgmilitant,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

######################################
## This section contains code for predicted probabilities.

## Territory

terr.min <- 0; terr.max <- 2

interact <- c(0, 1*mean(orgmilitant), 2*mean(orgmilitant))

logit.1.constant <- na.omit(data.frame(ORGST9, orgmilitant, domorgprot, orgpop, relorg, natorgnew, polity2, ORGST9*orgmilitant))

terr.iv <- matrix(apply(logit.1.constant,2,mean),nrow=3,ncol=8,byrow=T)
terr.iv[,1] <- seq(terr.min,terr.max,length.out=3)
terr.iv[,8] <- interact

## Simulation
set.seed(12591) ## For replicability purposes
S <- 1000
pprobs <- matrix(NA,S,3)
betas <- mvrnorm(S,coef(logit.1),VCV.logit.1)

for(s in 1:S)
{
  pprobs[s,] <- exp(cbind(1,terr.iv)%*%betas[s,])/(1+exp(cbind(1,terr.iv)%*%betas[s,]))
}

## CIs
cis.99 <- apply(pprobs,2,quantile,c(0.005,0.995))
cis.95 <- apply(pprobs,2,quantile,c(0.025,0.975))
cis.90 <- apply(pprobs,2,quantile,c(0.05,0.95))

colnames(terr.iv) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity", "ORGST9*orgmilitant")

terr.prob <- exp(cbind(1,terr.iv)%*%coef(logit.1))/(1+exp(cbind(1,terr.iv)%*%coef(logit.1)))

terr.plot <- qplot(xlab="", ylab="Predicted Probability of Negotiation") + 
geom_bar(aes(terr.iv[,1],terr.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(0:2), labels=c("No Territorial Control","Some Territorial Control","Full Territorial Control")) + 
scale_y_continuous(labels=percent, limits = c(0,.7)) + 
geom_segment(aes(x=0:2,xend=0:2,y=cis.95[1,],yend=cis.95[2,]),colour="black", size=2) + 
geom_segment(aes(x=0:2,xend=0:2,y=cis.90[1,],yend=cis.90[2,]),colour="black", size=3) + geom_segment(aes(x=0:2,xend=0:2,y=cis.99[1,],yend=cis.99[2,]),colour="black", size=1) + 
theme_bw() + 
theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold", hjust=.7), axis.title=element_text(size=28,face="bold"), axis.text.x = element_text(angle = 45, hjust=1))

## Territory & Fringe

terr.frin <- matrix(apply(logit.1.constant,2,mean),nrow=3,ncol=8,byrow=T)
terr.frin[,1] <- seq(terr.min,terr.max,length.out=3)
terr.frin[,4] <- 1
terr.frin[,8] <- interact

colnames(terr.frin) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity", "ORGST9*orgmilitant")

## Simulation
pprobs.frin <- matrix(NA,S,3)

for(s in 1:S)
{
  pprobs.frin[s,] <- exp(cbind(1,terr.frin)%*%betas[s,])/(1+exp(cbind(1,terr.frin)%*%betas[s,]))
}

## CIs
cis.99.frin <- apply(pprobs.frin,2,quantile,c(0.005,0.995))
cis.95.frin <- apply(pprobs.frin,2,quantile,c(0.025,0.975))
cis.90.frin <- apply(pprobs.frin,2,quantile,c(0.05,0.95))

terr.frin.prob <- exp(cbind(1,terr.frin)%*%coef(logit.1))/(1+exp(cbind(1,terr.frin)%*%coef(logit.1)))

terr.frin.plot <- qplot(xlab="", ylab="Predicted Probability of Negotiation") + geom_bar(aes(terr.frin[,1],terr.frin.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(0:2), labels=c("No Territorial Control","Some Territorial Control","Full Territorial Control")) + scale_y_continuous(labels=percent, limits = c(0,.9)) + geom_segment(aes(x=0:2,xend=0:2,y=cis.95.frin[1,],yend=cis.95.frin[2,]),colour="black", size=2) + geom_segment(aes(x=0:2,xend=0:2,y=cis.90.frin[1,],yend=cis.90.frin[2,]),colour="black", size=3) + geom_segment(aes(x=0:2,xend=0:2,y=cis.99.frin[1,],yend=cis.99.frin[2,]),colour="black", size=1) + theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold", hjust=.7), axis.title=element_text(size=28,face="bold"), axis.text.x = element_text(angle = 45, hjust=1))

## Territory & Competitive

terr.comp <- matrix(apply(logit.1.constant,2,mean),nrow=3,ncol=8,byrow=T)
terr.comp[,1] <- seq(terr.min,terr.max,length.out=3)
terr.comp[,4] <- 2
terr.comp[,8] <- interact

colnames(terr.comp) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity", "ORGST9*orgmilitant")

## Simulation
pprobs.comp <- matrix(NA,S,3)

for(s in 1:S)
{
  pprobs.comp[s,] <- exp(cbind(1,terr.comp)%*%betas[s,])/(1+exp(cbind(1,terr.comp)%*%betas[s,]))
}

## CIs
cis.99.comp <- apply(pprobs.comp,2,quantile,c(0.005,0.995))
cis.95.comp <- apply(pprobs.comp,2,quantile,c(0.025,0.975))
cis.90.comp <- apply(pprobs.comp,2,quantile,c(0.05,0.95))

terr.comp.prob <- exp(cbind(1,terr.comp)%*%coef(logit.1))/(1+exp(cbind(1,terr.comp)%*%coef(logit.1)))

terr.comp.plot <- qplot(xlab="", ylab="Predicted Probability of Negotiation") + geom_bar(aes(terr.comp[,1],terr.comp.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(0:2), labels=c("No Territorial Control","Some Territorial Control","Full Territorial Control")) + scale_y_continuous(labels=percent, limits = c(0,.9)) + geom_segment(aes(x=0:2,xend=0:2,y=cis.95.comp[1,],yend=cis.95.comp[2,]),colour="black", size=2) + geom_segment(aes(x=0:2,xend=0:2,y=cis.90.comp[1,],yend=cis.90.comp[2,]),colour="black", size=3) + geom_segment(aes(x=0:2,xend=0:2,y=cis.99.comp[1,],yend=cis.99.comp[2,]),colour="black", size=1) + theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold", hjust=.7), axis.title=element_text(size=28,face="bold"), axis.text.x = element_text(angle = 45, hjust=1))

## Territory & Dominant

terr.dom <- matrix(apply(logit.1.constant,2,mean),nrow=3,ncol=8,byrow=T)
terr.dom[,1] <- seq(terr.min,terr.max,length.out=3)
terr.dom[,4] <- 3
terr.dom[,8] <- interact

colnames(terr.dom) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity", "ORGST9*orgmilitant")

## Simulation
pprobs.dom <- matrix(NA,S,3)

for(s in 1:S)
{
  pprobs.dom[s,] <- exp(cbind(1,terr.dom)%*%betas[s,])/(1+exp(cbind(1,terr.dom)%*%betas[s,]))
}

## CIs
cis.99.dom <- apply(pprobs.dom,2,quantile,c(0.005,0.995))
cis.95.dom <- apply(pprobs.dom,2,quantile,c(0.025,0.975))
cis.90.dom <- apply(pprobs.dom,2,quantile,c(0.05,0.95))

terr.dom.prob <- exp(cbind(1,terr.dom)%*%coef(logit.1))/(1+exp(cbind(1,terr.dom)%*%coef(logit.1)))

terr.dom.plot <- qplot(xlab="", ylab="Predicted Probability of Negotiation") + geom_bar(aes(terr.dom[,1],terr.dom.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(0:2), labels=c("No Territorial Control","Some Territorial Control","Full Territorial Control")) + scale_y_continuous(labels=percent, limits = c(0,.9)) + geom_segment(aes(x=0:2,xend=0:2,y=cis.95.dom[1,],yend=cis.95.dom[2,]),colour="black", size=2) + geom_segment(aes(x=0:2,xend=0:2,y=cis.90.dom[1,],yend=cis.90.dom[2,]),colour="black", size=3) + geom_segment(aes(x=0:2,xend=0:2,y=cis.99.dom[1,],yend=cis.99.dom[2,]),colour="black", size=1) + theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold", hjust=.7), axis.title=element_text(size=28,face="bold"), axis.text.x = element_text(angle = 45, hjust=1))

## Militancy

mil.min <- 0 ; mil.max <- 1

interact2 <- c(0,mean(ORGST9))

mil.iv <- matrix(apply(logit.1.constant,2,mean),nrow=2,ncol=8,byrow=T)
mil.iv[,2] <- seq(mil.min,mil.max,length.out=2)
mil.iv[,8] <- interact2

pprobs.mil <- matrix(NA,S,2)
betas.mil <- mvrnorm(S,coef(logit.1),VCV.logit.1)

for(s in 1:S)
{
  pprobs.mil[s,] <- exp(cbind(1,mil.iv)%*%betas.mil[s,])/(1+exp(cbind(1,mil.iv)%*%betas.mil[s,]))
}

## CIs
cis.99.mil <- apply(pprobs.mil,2,quantile,c(0.005,0.995))
cis.95.mil <- apply(pprobs.mil,2,quantile,c(0.025,0.975))
cis.90.mil <- apply(pprobs.mil,2,quantile,c(0.05,0.95))

colnames(mil.iv) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

mil.prob <- exp(cbind(1,mil.iv)%*%coef(logit.1))/(1+exp(cbind(1,mil.iv)%*%coef(logit.1)))

mil.plot <- qplot(xlab="", ylab="Predicted Probability of Negotiation") + geom_bar(aes(mil.iv[,2],mil.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(0:1), labels=c("No Militant Wing","Militant Wing")) + scale_y_continuous(labels=percent, limits = c(0,.7)) + geom_segment(aes(x=0:1,xend=0:1,y=cis.95.mil[1,],yend=cis.95.mil[2,]),colour="black", size=2) + geom_segment(aes(x=0:1,xend=0:1,y=cis.90.mil[1,],yend=cis.90.mil[2,]),colour="black", size=3) + geom_segment(aes(x=0:1,xend=0:1,y=cis.99.mil[1,],yend=cis.99.mil[2,]),colour="black", size=1) + theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold"), axis.title=element_text(size=28,face="bold"))

## Protest

prot.min <- 0 ; prot.max <- 5

interact3 <- mean(ORGST9) * mean(orgmilitant)

prot.iv <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)
prot.iv[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv[,8] <- interact3

pprobs.prot <- matrix(NA,S,6)
betas.prot <- mvrnorm(S,coef(logit.1),VCV.logit.1)

for(s in 1:S)
{
  pprobs.prot[s,] <- exp(cbind(1,prot.iv)%*%betas.prot[s,])/(1+exp(cbind(1,prot.iv)%*%betas.prot[s,]))
}

## CIs
cis.99.prot <- apply(pprobs.prot,2,quantile,c(0.005,0.995))
cis.95.prot <- apply(pprobs.prot,2,quantile,c(0.025,0.975))
cis.90.prot <- apply(pprobs.prot,2,quantile,c(0.05,0.95))

colnames(prot.iv) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2","ORGST9*orgmilitant")

prot.prob <- exp(cbind(1,prot.iv)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv)%*%coef(logit.1)))

prot.plot <- qplot(xlab="", ylab="Predicted Probability of Negotiation") + geom_bar(aes(prot.iv[,3],prot.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(0:5), labels=c("No Opposition", "Verbal Opposition", "Symbolic Resistance", "Small Demonstrations", "Medium Demonstrations", "Large Demonstrations")) + theme(axis.text.x=element_text(angle = -30, hjust = 0)) + scale_y_continuous(labels=percent, limits = c(0,.7)) + geom_segment(aes(x=0:5,xend=0:5,y=cis.95.prot[1,],yend=cis.95.prot[2,]),colour="black", size=2) + geom_segment(aes(x=0:5,xend=0:5,y=cis.90.prot[1,],yend=cis.90.prot[2,]),colour="black", size=3) + geom_segment(aes(x=0:5,xend=0:5,y=cis.99.prot[1,],yend=cis.99.prot[2,]),colour="black", size=1) + theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold", hjust=.7), axis.title=element_text(size=28,face="bold"),axis.text.x = element_text(angle = 45, hjust=1))

## Popularity

pop.min <- 1 ; pop.max <- 3

pop.iv <- matrix(apply(logit.1.constant,2,mean),nrow=3,ncol=8,byrow=T)
pop.iv[,1] <- 0
pop.iv[,4] <- seq(pop.min,pop.max,length.out=3)
pop.iv[,8] <- 0

pprobs.pop <- matrix(NA,S,3)
betas.pop <- mvrnorm(S,coef(logit.1),VCV.logit.1)

for(s in 1:S)
{
  pprobs.pop[s,] <- exp(cbind(1,pop.iv)%*%betas.pop[s,])/(1+exp(cbind(1,pop.iv)%*%betas.pop[s,]))
}

## CIs
cis.99.pop <- apply(pprobs.pop,2,quantile,c(0.005,0.995))
cis.95.pop <- apply(pprobs.pop,2,quantile,c(0.025,0.975))
cis.90.pop <- apply(pprobs.pop,2,quantile,c(0.05,0.95))

colnames(pop.iv) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

pop.prob <- exp(cbind(1,pop.iv)%*%coef(logit.1))/(1+exp(cbind(1,pop.iv)%*%coef(logit.1)))

pop.plot <- qplot(xlab="",ylab="Predicted Probability of Negotiation") + geom_bar(aes(pop.iv[,4],pop.prob), stat = "identity", fill = "darkgrey") + scale_x_continuous(breaks=c(1:3), labels=c("No Support", "Competing Group", "Dominant Group")) + scale_y_continuous(labels=percent, limits = c(0,.7)) + geom_segment(aes(x=1:3,xend=1:3,y=cis.95.pop[1,],yend=cis.95.pop[2,]),colour="black", size=2) + geom_segment(aes(x=1:3,xend=1:3,y=cis.90.pop[1,],yend=cis.90.pop[2,]),colour="black", size=3) + geom_segment(aes(x=1:3,xend=1:3,y=cis.99.pop[1,],yend=cis.99.pop[2,]),colour="black", size=1) + theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text=element_text(size=18, face="bold"), axis.title=element_text(size=28,face="bold"))

######################################
## This section contains the code for an autoregressive model.

## Lag function
lg <- function(x)c(NA, x[1:(length(x)-1)])

MAROB.Negot.lag <- ddply(MAROB.Negot, ~orgId, transform, ACTIVENEGOTIATION_lag = lg(ACTIVENEGOTIATION))

MAROB.Negot.lag <- na.omit(MAROB.Negot.lag)

logit.2 <- glm(ACTIVENEGOTIATION ~ ORGST9 + orgmilitant + domorgprot + orgpop + relorg + natorgnew + polity2 + ORGST9*orgmilitant + ACTIVENEGOTIATION_lag, data = MAROB.Negot.lag, family = "binomial")

summary(logit.2)

## Clustered SE
cluster.vcv.2 <- cluster.vcov(logit.2, cluster = MAROB.Negot.lag$orgId)

## Results
new.coef.2 <- coeftest(logit.2, cluster.vcv.2)

## Check the interaction term
VCV.logit.2 <- cluster.vcv.2

P2 <- length(coef(logit.2))

CI <- .95

s.orgmilitant <- seq(min(orgmilitant),max(orgmilitant),length.out=1041)

## The marginal effect of territorial control
Marg.ORGST9.2 <- coef(logit.2)[2] + coef(logit.2)[10]*s.orgmilitant

SE.Marg.ORGST9.2 <- sqrt(VCV.logit.2[2,2] + (s.orgmilitant^2)*VCV.logit.2[10,10] + 2*s.orgmilitant*VCV.logit.2[2,10])

Lowerci.ORGST9.2 <- Marg.ORGST9.2 - qt((1-CI)/2,1041-P)*SE.Marg.ORGST9.2
Upperci.ORGST9.2 <- Marg.ORGST9.2 + qt((1-CI)/2,1041-P)*SE.Marg.ORGST9.2

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

plot(s.orgmilitant,Marg.ORGST9.2,type="n",
     xlab="Observed Values of Militancy",
     ylab="Marginal Effect of Territorial Control",
     ylim=c(min(Lowerci.ORGST9.2)-2,max(Upperci.ORGST9.2)+2))

lines(s.orgmilitant,Marg.ORGST9.2,lwd=5,lty=1)
lines(s.orgmilitant,Lowerci.ORGST9.2,lwd=4,lty=2, col="blue")
lines(s.orgmilitant,Upperci.ORGST9.2,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

############################
## Model accounting for lagged territorial control

MAROB.Negot.lag2 <- ddply(MAROB.Negot, ~orgId, transform, ORGST9_lag = lg(ORGST9))

MAROB.Negot.lag2 <- na.omit(MAROB.Negot.lag2)

logit.3 <- glm(ACTIVENEGOTIATION ~ ORGST9_lag + orgmilitant + domorgprot + orgpop + relorg + natorgnew + polity2+ ORGST9_lag*orgmilitant, data = MAROB.Negot.lag2, family = "binomial")

## Clustered SE
cluster.vcv.3 <- cluster.vcov(logit.3, cluster = MAROB.Negot.lag2$orgId)

## Results
new.coef.3 <- coeftest(logit.3, cluster.vcv.3)

summary(logit.3)

## Check the interaction term
VCV.logit.3 <- cluster.vcv.3

P3 <- length(coef(logit.3))

CI <- .95

s.orgmilitant <- seq(min(orgmilitant),max(orgmilitant),length.out=1041)

## The marginal effect of territorial control
Marg.ORGST9.3 <- coef(logit.3)[2] + coef(logit.3)[9]*s.orgmilitant

SE.Marg.ORGST9.3 <- sqrt(VCV.logit.3[2,2] + (s.orgmilitant^2)*VCV.logit.3[9,9] + 2*s.orgmilitant*VCV.logit.3[2,9])

Lowerci.ORGST9.3 <- Marg.ORGST9.3 - qt((1-CI)/2,1041-P)*SE.Marg.ORGST9.3
Upperci.ORGST9.3 <- Marg.ORGST9.3 + qt((1-CI)/2,1041-P)*SE.Marg.ORGST9.3

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

plot(s.orgmilitant,Marg.ORGST9.3,type="n",
     xlab="Observed Values of Militancy",
     ylab="Marginal Effect of Territorial Control",
     ylim=c(min(Lowerci.ORGST9.3)-2,max(Upperci.ORGST9.3)+2))

lines(s.orgmilitant,Marg.ORGST9.3,lwd=5,lty=1)
lines(s.orgmilitant,Lowerci.ORGST9.3,lwd=4,lty=2, col="blue")
lines(s.orgmilitant,Upperci.ORGST9.3,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

################################
## This is a test when we collapse territorial control to become binary

## Create the binary measure
MAROB.Negot$ORGST9_b <- ifelse(ORGST9 >=1,1,0)

## Run the test
logit.4 <- glm(ACTIVENEGOTIATION ~ ORGST9_b + orgmilitant + domorgprot + orgpop + relorg + natorgnew + polity2 + ORGST9_b*orgmilitant, data = MAROB.Negot, family = "binomial")

## Clustered SE
cluster.vcv.4 <- cluster.vcov(logit.4, cluster = orgId)

## Results
new.coef.4 <- coeftest(logit.4, cluster.vcv.4)

summary(logit.4)

## Check the interaction term
VCV.logit.4 <- cluster.vcv.4

P4 <- length(coef(logit.4))

CI <- .95

s.orgmilitant <- seq(min(orgmilitant),max(orgmilitant),length.out=1135)

## The marginal effect of territorial control
Marg.ORGST9.4 <- coef(logit.4)[2] + coef(logit.4)[9]*s.orgmilitant

SE.Marg.ORGST9.4 <- sqrt(VCV.logit.4[2,2] + (s.orgmilitant^2)*VCV.logit.4[9,9] + 2*s.orgmilitant*VCV.logit.4[2,9])

Lowerci.ORGST9.4 <- Marg.ORGST9.4 - qt((1-CI)/2,1135-P)*SE.Marg.ORGST9.4
Upperci.ORGST9.4 <- Marg.ORGST9.4 + qt((1-CI)/2,1135-P)*SE.Marg.ORGST9.4

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

plot(s.orgmilitant,Marg.ORGST9.4,type="n",
     xlab="Observed Values of Militancy",
     ylab="Marginal Effect of Territorial Control",
     ylim=c(min(Lowerci.ORGST9.4)-2,max(Upperci.ORGST9.4)+2))

lines(s.orgmilitant,Marg.ORGST9.4,lwd=5,lty=1)
lines(s.orgmilitant,Lowerci.ORGST9.4,lwd=4,lty=2, col="blue")
lines(s.orgmilitant,Upperci.ORGST9.4,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

##################################
## This is a test without ideological controls.

logit.5 <- glm(ACTIVENEGOTIATION ~ ORGST9 + orgmilitant + domorgprot + orgpop+ polity2 + ORGST9*orgmilitant, data = MAROB.Negot, family = "binomial")

## Clustered SE
cluster.vcv.5 <- cluster.vcov(logit.5, cluster = orgId)

## Results
new.coef.5 <- coeftest(logit.5, cluster.vcv.5)

summary(logit.5)

## Check the interaction term
VCV.logit.5 <- cluster.vcv.5

P5 <- length(coef(logit.5))

CI <- .95

s.orgmilitant <- seq(min(orgmilitant),max(orgmilitant),length.out=1135)

## The marginal effect of territorial control
Marg.ORGST9.5 <- coef(logit.5)[2] + coef(logit.5)[7]*s.orgmilitant

SE.Marg.ORGST9.5 <- sqrt(VCV.logit.5[2,2] + (s.orgmilitant^2)*VCV.logit.5[7,7] + 2*s.orgmilitant*VCV.logit.5[2,7])

Lowerci.ORGST9.5 <- Marg.ORGST9.5 - qt((1-CI)/2,1135-P)*SE.Marg.ORGST9.5
Upperci.ORGST9.5 <- Marg.ORGST9.5 + qt((1-CI)/2,1135-P)*SE.Marg.ORGST9.5

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

plot(s.orgmilitant,Marg.ORGST9.5,type="n",
     xlab="Observed Values of Militancy",
     ylab="Marginal Effect of Territorial Control",
     ylim=c(min(Lowerci.ORGST9.5)-2,max(Upperci.ORGST9.5)+2))

lines(s.orgmilitant,Marg.ORGST9.5,lwd=5,lty=1)
lines(s.orgmilitant,Lowerci.ORGST9.5,lwd=4,lty=2, col="blue")
lines(s.orgmilitant,Upperci.ORGST9.5,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

#########################
## This section substitutes terrorist tactics for militancy to better engage with Thomas (2014)

logit.6 <- glm(ACTIVENEGOTIATION ~ ORGST9 + ORGST7 + domorgprot + orgpop + relorg + natorgnew + polity2 + ORGST9*ORGST7, data = MAROB.Negot, family = "binomial")

## Clustered SE
cluster.vcv.6 <- cluster.vcov(logit.6, cluster = orgId)

## Results
new.coef.6 <- coeftest(logit.6, cluster.vcv.6)

summary(logit.6)

## Check the interaction term
VCV.logit.6 <- cluster.vcv.6

P6 <- length(coef(logit.6))

CI <- .95

s.ORGST7 <- seq(min(ORGST7),max(ORGST7),length.out=1135)

## The marginal effect of territorial control
Marg.ORGST9.6 <- coef(logit.6)[2] + coef(logit.6)[9]*s.ORGST7

SE.Marg.ORGST9.6 <- sqrt(VCV.logit.6[2,2] + (s.ORGST7^2)*VCV.logit.6[9,9] + 2*s.ORGST7*VCV.logit.6[2,9])

Lowerci.ORGST9.6 <- Marg.ORGST9.6 - qt((1-CI)/2,1135-P)*SE.Marg.ORGST9.6
Upperci.ORGST9.6 <- Marg.ORGST9.6 + qt((1-CI)/2,1135-P)*SE.Marg.ORGST9.6

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

plot(s.ORGST7,Marg.ORGST9.6,type="n",
     xlab="Observed Values of Terrorism",
     ylab="Marginal Effect of Territorial Control",
     ylim=c(min(Lowerci.ORGST9.6)-2,max(Upperci.ORGST9.6)+2))

lines(s.ORGST7,Marg.ORGST9.6,lwd=5,lty=1)
lines(s.ORGST7,Lowerci.ORGST9.6,lwd=4,lty=2, col="blue")
lines(s.ORGST7,Upperci.ORGST9.6,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

####
## Checking the marginal effect of terrorism

s.ORGST9 <- seq(min(ORGST9),max(ORGST9),length.out=1135)

## The marginal effect of terrorism for varying values of territorial control.
Marg.ORGST7 <- coef(logit.6)[3] + coef(logit.6)[9]*s.ORGST9

SE.Marg.ORGST7 <- sqrt(VCV.logit.6[3,3] + (s.ORGST9^2)*VCV.logit.6[9,9] + 2*s.ORGST9*VCV.logit.6[3,9])

Lowerci.ORGST7 <- Marg.ORGST7 - qt((1-CI)/2,1135-P)*SE.Marg.ORGST7
Upperci.ORGST7 <- Marg.ORGST7 + qt((1-CI)/2,1135-P)*SE.Marg.ORGST7

## Plot for visualization. N.B. the only points that exist on this plot occur at territorial control = 0,1,2. This figure will not be presented.
plot(s.ORGST9,Marg.ORGST7,type="n",
     xlab="Observed Values of Territorial Control",
     ylab="Marginal Effect of Terrorism",
     ylim=c(min(Lowerci.ORGST7)-2,max(Upperci.ORGST7)+2))

lines(s.ORGST9,Marg.ORGST7,lwd=5,lty=1)
lines(s.ORGST9,Lowerci.ORGST7,lwd=4,lty=2, col="blue")
lines(s.ORGST9,Upperci.ORGST7,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

#########################
## Country and year fixed effects

logit.fe <- glm(ACTIVENEGOTIATION ~ ORGST9 + orgmilitant + domorgprot + orgpop + relorg + natorgnew + polity2 + ORGST9*orgmilitant + as.factor(country) + as.factor(year) , data = MAROB.Negot, family = "binomial")

## Clustered SE
cluster.vcv.fe <- cluster.vcov(logit.fe, cluster = orgId)

## Results
new.coef.fe <- coeftest(logit.fe, cluster.vcv.fe)

summary(logit.fe)

## Check the interaction term
VCV.logit.fe <- cluster.vcv.fe

Pfe <- length(coef(logit.fe))

CI <- .95

s.orgmilitant <- seq(min(orgmilitant),max(orgmilitant),length.out=1135)

## The marginal effect of territorial control
Marg.ORGST9.fe <- coef(logit.fe)[2] + coef(logit.fe)[44]*s.orgmilitant

SE.Marg.ORGST9.fe <- sqrt(VCV.logit.fe[2,2] + (s.orgmilitant^2)*VCV.logit.fe[44,44] + 2*s.orgmilitant*VCV.logit.fe[2,9])

Lowerci.ORGST9.fe <- Marg.ORGST9.fe - qt((1-CI)/2,1135-Pfe)*SE.Marg.ORGST9.fe
Upperci.ORGST9.fe <- Marg.ORGST9.fe + qt((1-CI)/2,1135-Pfe)*SE.Marg.ORGST9.fe

## Plot for visualization. N.B. that the only true points on this plot occur at orgmilitant=0 and orgmilitant=1. The rest is only for visualization and will not be presented.

plot(s.orgmilitant,Marg.ORGST9.fe,type="n",
     xlab="Observed Values of Militancy",
     ylab="Marginal Effect of Territorial Control",
     ylim=c(min(Lowerci.ORGST9.fe)-3,max(Upperci.ORGST9.fe)+2))

lines(s.orgmilitant,Marg.ORGST9.fe,lwd=5,lty=1)
lines(s.orgmilitant,Lowerci.ORGST9.fe,lwd=4,lty=2, col="blue")
lines(s.orgmilitant,Upperci.ORGST9.fe,lwd=4,lty=2, col="blue")

abline(h=0,lwd=2)

#########################
## This produces a table of descriptive statistics for variables used in the analysis
vars <- cbind(ACTIVENEGOTIATION, ORGST9, orgmilitant, ORGST7, domorgprot, orgpop, relorg, natorgnew, polity2)

options(scipen=100)
options(digits=2)
desc.table <- stat.desc(vars)

## Transpose table
desc.table <- t(desc.table)
xtable(desc.table)

########################
## This code will create our main table of results

texreg(list(new.coef, new.coef.4, new.coef.6, new.coef.3, new.coef.2, new.coef.fe), dcolumn=T, stars=.05)

#########################
## This section will look at the combined effects of territorial control, organizational popularity, and domestic protest

prot.iv.c1 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c1) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c1[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c1[,1] <- 0
prot.iv.c1[,4] <- 1
prot.iv.c1[,8] <- interact3

prot.prob.c1 <- exp(cbind(1,prot.iv.c1)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c1)%*%coef(logit.1)))

##

prot.iv.c2 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c2) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c2[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c2[,1] <- 0
prot.iv.c2[,4] <- 2
prot.iv.c2[,8] <- interact3

prot.prob.c2 <- exp(cbind(1,prot.iv.c2)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c2)%*%coef(logit.1)))

##

prot.iv.c3 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c3) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c3[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c3[,1] <- 0
prot.iv.c3[,4] <- 3
prot.iv.c3[,8] <- interact3

prot.prob.c3 <- exp(cbind(1,prot.iv.c3)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c3)%*%coef(logit.1)))

##

prot.iv.c4 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c4) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c4[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c4[,1] <- 1
prot.iv.c4[,4] <- 1
prot.iv.c4[,8] <- interact3

prot.prob.c4 <- exp(cbind(1,prot.iv.c4)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c4)%*%coef(logit.1)))

##

prot.iv.c5 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c5) <-c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c5[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c5[,1] <- 1
prot.iv.c5[,4] <- 2
prot.iv.c5[,8] <- interact3

prot.prob.c5 <- exp(cbind(1,prot.iv.c5)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c5)%*%coef(logit.1)))

##

prot.iv.c6 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c6) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c6[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c6[,1] <- 1
prot.iv.c6[,4] <- 3
prot.iv.c6[,8] <- interact3

prot.prob.c6 <- exp(cbind(1,prot.iv.c6)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c6)%*%coef(logit.1)))

##

prot.iv.c7 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c7) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c7[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c7[,1] <- 2
prot.iv.c7[,4] <- 1
prot.iv.c7[,8] <- interact3

prot.prob.c7 <- exp(cbind(1,prot.iv.c7)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c7)%*%coef(logit.1)))

##

prot.iv.c8 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c8) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c8[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c8[,1] <- 2
prot.iv.c8[,4] <- 2
prot.iv.c8[,8] <- interact3

prot.prob.c8 <- exp(cbind(1,prot.iv.c8)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c8)%*%coef(logit.1)))

##

prot.iv.c9 <- matrix(apply(logit.1.constant,2,mean),nrow=6,ncol=8,byrow=T)

colnames(prot.iv.c9) <- c("ORGST9", "orgmilitant", "domorgprot", "orgpop", "relorg", "natorgnew", "polity2", "ORGST9*orgmilitant")

prot.iv.c9[,3] <- seq(prot.min,prot.max,length.out=6)
prot.iv.c9[,1] <- 2
prot.iv.c9[,4] <- 3
prot.iv.c9[,8] <- interact3

prot.prob.c9 <- exp(cbind(1,prot.iv.c9)%*%coef(logit.1))/(1+exp(cbind(1,prot.iv.c9)%*%coef(logit.1)))

##

combined.pp <- cbind(prot.prob.c1,prot.prob.c2,prot.prob.c3,prot.prob.c4,prot.prob.c5,prot.prob.c6,prot.prob.c7,prot.prob.c8,prot.prob.c9)

combined.pp.shape = melt(combined.pp)

combined_heatmap <- ggplot(combined.pp.shape, aes(X2,X1))+
scale_fill_gradient("Predicted\nProbability",low="white",high="black", limits=c(0,1), breaks=c(0,0.5,1))+
geom_tile(aes(fill=value))+
scale_x_continuous(breaks=c(1:9), labels = c(1,2,3,1,2,3,1,2,3)) + 
scale_y_reverse(breaks=c(1:6), labels=0:5) + 
xlab("Organizational Popularity") + 
ylab("Domestic Protest") + 
annotate("text", label = "0", x=2, y=0,size=8, fontface="bold") + 
annotate("text", label = "1", x=5, y=0, size=8, fontface="bold") + 
annotate("text", label = "2", x=8, y=0, size=8, fontface="bold") + 
annotate("text", label = "Territorial Control", x=5, y = -.5, size = 8, fontface="bold") + 
theme_bw() + 
theme(legend.title=element_text(size=18,face="bold"), legend.text=element_text(size=18, face="bold"), legend.key.size = unit(2.5, "cm"), legend.position="bottom", panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), panel.grid.major.x=element_blank(), panel.grid.major.y=element_blank(), axis.text=element_text(size=18, face="bold", hjust=.7), axis.title=element_text(size=28,face="bold")) +
geom_segment(aes(x=.5,xend=.5,y=.5,yend=6.5)) + 
geom_segment(aes(x=3.5,xend=3.5,y=.5,yend=6.5)) + 
geom_segment(aes(x=6.5,xend=6.5,y=.5,yend=6.5)) + 
geom_segment(aes(x=.5,xend=9.5,y=6.5,yend=6.5)) + 
geom_segment(aes(x=.5,xend=9.5,y=.5,yend=.5)) + 
geom_segment(aes(x=9.5,xend=9.5,y=6.5,yend=.5))

########################